setClass("StdbwIClusterScore",contains="InternalClusterScore",
		prototype=prototype(
				.description="Stdbw Cluster Internal Score class"
		
		)
)  

#global functions

#helper functions

calculateVarVec <- function(points){
	variance <- ( unlist( lapply(1:dim(points)[2] ,FUN=function(n,set){
									#this is slow var(set[,n])
									vec <-set[,n]
									len <- length(vec)
									if(len==1){return(0)}
									{sum(vec*vec) - sum(vec)*sum(vec)/len}/{len-1}
									#var(set[,n])
								},
								points ) ) )
	return(variance)
}

#function to calculate density around given center point
#INPUT:
#	center -center point of cluster
#	Rn - radius of neighbourhood
#	points - points

calDens <- function(center,Rn,points){
	exPoints <- rbind(center,points)
	#distances to given center
	distances <- as.matrix(dist(exPoints))[1]
	
	numOfNeigh <- length(which(distances<=Rn))
	if(numOfNeigh==0){ warning("No neighbours found");return(as.integer(1))}
	return(numOfNeigh)
}
calDens<-cmpfun(calDens)

#methods

#calclustes silhouette index for given clusterization
setMethod("ScoreSet",
		signature="StdbwIClusterScore",
		definition=function(.Object,inputSet,clusters,...){
			#check for X and clusters
			callNextMethod(.Object,inputSet,clusters,...)
			if(isOneCluster(clusters))return(-.Machine$double.xmax)
			
			#calclulate number of clusters
			c <- length(unique(clusters))
			
			#Calulate vector of variances in Data set
			dataVarianceVec <- calculateVarVec(inputSet$X)
			dataVarianceVecLen <- sqrt( sum(dataVarianceVec^2 ) )
			
			#stopifnot(!is.nan(dataVarianceVec))
			#stopifnot(!is.nan(dataVarianceVecLen))
			
			#Calculate vectors of variances in clusters
			#each row contains variance vector
			#print("cluster Variance")
			clustersVarianceMatrix <- matrix( unlist( lapply(1:c,FUN=function(n,points,clusters){
											clusterElements <- which(clusters==n)
											varVec <- calculateVarVec(points[clusterElements,,drop=FALSE])							
										} ,inputSet$X,clusters ) )
									, nrow=c,byrow=TRUE)
							
			#stopifnot(!is.nan(clustersVarianceMatrix))
			
			#calculate lengths fo clusters variance vectors
			#each field of vector contains length of variance vector
			clustersVarianceLen <-unlist(lapply(1:c,FUN=function(n){
								sqrt( sum( clustersVarianceMatrix[n,,drop=FALSE]*clustersVarianceMatrix[n,,drop=FALSE] ) )
							},...) )
			
			#stopifnot(!is.nan(clustersVarianceLen))
			
			#calculate cluster centers
			#each row contains cluster center
			cluCentsMatrix <- matrix(unlist(lapply(1:c,FUN=function(n,clusters,points){										
										clusterElements <- which(clusters==n)
										clusterElementsNumber <- length(clusterElements)
										clusterCenter <- colSums(points[clusterElements,,drop=FALSE])/clusterElementsNumber							
									},clusters,inputSet$X)
					),nrow=c,byrow=TRUE)
			
			#stopifnot(!is.nan(cluCentsMatrix))

			#calculate length of center vector of cluster
#			clustersCenterLengths <-unlist(lapply(1:c,FUN=function(n){
#								sqrt( sum( (cluCentsMatrix[n,])^2 ) )
#							},...) )
			
			#stopifnot(!is.nan(clustersCenterLengths))

			#calculate scatter
			scatter <-  sum(clustersVarianceLen)/{c*dataVarianceVecLen}
			
			#stopifnot(!is.nan(scatter))
			
			#calculate radius of neighbourhood
			Rnei <- sqrt( sum( clustersVarianceLen ) )/c
			
			
			#stopifnot(!is.nan(Rnei))

			#make centers grid
			gri <- expand.grid(1:c,1:c)
			#upper triangle of grid matrix
			gri <- gri[which(gri$Var1>gri$Var2),]
			
			#calculate densSum
			densSum <- sum(apply(gri,1,FUN=function(gp){								
								numerator <- calDens( {cluCentsMatrix[gp[1],] +  cluCentsMatrix[gp[2],]}/2 ,Rnei,inputSet$X)
								denominator <-max(
												calDens(  cluCentsMatrix[gp[1],]  ,Rnei ,inputSet$X) ,
												calDens( cluCentsMatrix[gp[2],] ,Rnei,inputSet$X)
												)
								return(numerator/denominator)
							}))
			
			#stopifnot(!is.nan(densSum))
			#calculate densBW
			densBW <- densSum/{c*{c-1}}
			
			#stopifnot(!is.nan(densBW))
			#calculate stdbw
			stdbw <- scatter + densBW
			
			#stopifnot(!is.nan(stdbw))
			#print("STDBW: ");print(-stdbw)
			return(-stdbw) #maximize
		}
)  

